Chapter 6 Does the antibiotic administration change the microbial community?

6.1 Alpha diversity

alpha_div_meta <- alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == sample))%>%
  filter(time_point %in% c("Acclimation", "Antibiotics") ) 

label_map <- c(
  "richness" = "Species Richness",
  "neutral" = "Neutral Diversity",
  "phylogenetic" = "Phylogenetic Diversity")

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = "sample") %>%
  mutate(metric = factor(metric, levels = c("richness", "neutral", "phylogenetic")),
  time_point = factor(time_point, levels = c("Wild", "Acclimation"))) %>%
  filter(time_point %in% c("Wild", "Acclimation")) %>% 
  ggplot(aes(y = value, x = time_point, color=Population, fill=Population)) +
  geom_boxplot(width = 0.5,alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
  geom_jitter(width = 0.2,show.legend = FALSE) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#008080', "#d57d2c")) +
  facet_grid(metric ~ Population, scales = "free_y", labeller = labeller(metric = label_map, type = label_map))+
  theme_classic() +
  theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.text = element_text(size = 14),
    axis.title = element_text(size = 14),
    strip.text = element_text(size = 18, lineheight = 0.6),
      ) +
  ylab("Alpha diversity")

Richness

Modelq0GLMMNB <- glmer.nb(richness ~ time_point*Population+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(5.1645)  ( log )
Formula: richness ~ time_point * Population + (1 | individual)
   Data: alpha_div_meta

     AIC      BIC   logLik deviance df.resid 
   434.8    446.1   -211.4    422.8       43 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.61979 -0.56608  0.06074  0.51882  2.39923 

Random effects:
 Groups     Name        Variance Std.Dev.
 individual (Intercept) 0.2193   0.4683  
Number of obs: 49, groups:  individual, 27

Fixed effects:
                                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            3.7459     0.1623  23.073  < 2e-16 ***
time_pointAntibiotics                 -1.1523     0.1862  -6.190 6.02e-10 ***
PopulationWarm                         0.7055     0.2719   2.595  0.00946 ** 
time_pointAntibiotics:PopulationWarm  -0.3586     0.3037  -1.181  0.23769    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) tm_pnA PpltnW
tm_pntAntbt -0.466              
PopulatnWrm -0.598  0.278       
tm_pntAn:PW  0.292 -0.611 -0.461

Neutral

Model_neutral <- nlme::lme(fixed = neutral ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_neutral)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC    logLik
  342.5018 353.3418 -165.2509

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:    3.282619 7.923552

Fixed effects:  neutral ~ time_point * Population 
                                         Value Std.Error DF   t-value p-value
(Intercept)                           19.11713  2.078908 25  9.195758   0e+00
time_pointAntibiotics                -11.46317  2.832826 20 -4.046551   6e-04
PopulationWarm                        24.97345  3.534826 25  7.064972   0e+00
time_pointAntibiotics:PopulationWarm -20.91701  4.793363 20 -4.363745   3e-04
 Correlation: 
                                     (Intr) tm_pnA PpltnW
time_pointAntibiotics                -0.633              
PopulationWarm                       -0.588  0.373       
time_pointAntibiotics:PopulationWarm  0.374 -0.591 -0.632

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.8258256 -0.5685215 -0.1867832  0.5823712  2.0651455 

Number of Observations: 49
Number of Groups: 27 

Phylogenetic

Model_phylo <- nlme::lme(fixed = phylogenetic ~ time_point*Population, data = alpha_div_meta,
               random = ~ 1 | individual)
summary(Model_phylo)

Linear mixed-effects model fit by REML
  Data: alpha_div_meta 
       AIC      BIC   logLik
  203.6546 214.4946 -95.8273

Random effects:
 Formula: ~1 | individual
        (Intercept) Residual
StdDev:   0.7173285 1.688356

Fixed effects:  phylogenetic ~ time_point * Population 
                                         Value Std.Error DF   t-value p-value
(Intercept)                           5.365728 0.4446271 25 12.067930  0.0000
time_pointAntibiotics                -0.761455 0.6038652 20 -1.260969  0.2218
PopulationWarm                        1.097291 0.7560383 25  1.451370  0.1591
time_pointAntibiotics:PopulationWarm -0.702255 1.0216421 20 -0.687379  0.4997
 Correlation: 
                                     (Intr) tm_pnA PpltnW
time_pointAntibiotics                -0.631              
PopulationWarm                       -0.588  0.371       
time_pointAntibiotics:PopulationWarm  0.373 -0.591 -0.629

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.85993684 -0.51436620 -0.05130501  0.48944142  1.80209573 

Number of Observations: 49
Number of Groups: 27 

6.2 Beta diversity

samples_to_keep <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics")) %>% 
  select(Tube_code) %>% 
  pull()
subset_meta <- sample_metadata %>%
  filter(time_point %in% c("Acclimation", "Antibiotics"))

Richness

richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
               colnames(richness) %in% samples_to_keep])

betadisper(richness, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.03530 0.035300 10.073    999  0.002 **
Residuals 47 0.16471 0.003505                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(richness ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	1.9348458	0.1002347	6.065418	0.001
Population	1	2.1356198	0.1106358	6.694811	0.001
time_point:Population	1	0.8778553	0.0454773	2.751930	0.004
Residual	45	14.3548318	0.7436522	NA	NA
Total	48	19.3031527	1.0000000	NA	NA

Neutral

neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
               colnames(neutral) %in% samples_to_keep])

betadisper(neutral, subset_meta$time_point) %>% permutest(.) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.051657 0.051657 9.9224    999  0.005 **
Residuals 47 0.244689 0.005206                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(neutral ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	2.0429592	0.11061666	7.254623	0.001
Population	1	2.8764490	0.15574623	10.214376	0.001
time_point:Population	1	0.8770558	0.04748846	3.114457	0.003
Residual	45	12.6723552	0.68614864	NA	NA
Total	48	18.4688192	1.00000000	NA	NA

Phylogenetic

phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
               colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
Groups     1 0.67439 0.67439 55.932    999  0.001 ***
Residuals 47 0.56670 0.01206                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

adonis2(phylo ~ time_point*Population,
        data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
        permutations = 999,
        strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual),
        by="terms") %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
time_point	1	1.8783556	0.23032747	16.358708	0.001
Population	1	0.7551942	0.09260332	6.577030	0.001
time_point:Population	1	0.3545686	0.04347786	3.087959	0.023
Residual	45	5.1670341	0.63359135	NA	NA
Total	48	8.1551525	1.00000000	NA	NA

dbRDA

#Richness
cca_ord <- capscale(formula = richness ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Neutral
cca_ord <- capscale(formula = neutral ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(Population, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
    scale_fill_manual(values=c('#00808050', "#d57d2c50"))+
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

#Phylogenetic
cca_ord <- capscale(formula = phylo ~ subset_meta$time_point* subset_meta$Population)
CAP_df <- as.data.frame(vegan::scores(cca_ord, display = "sites")) %>%
  rownames_to_column('Tube_code') %>%
  left_join(subset_meta, by = 'Tube_code') %>%
  column_to_rownames('Tube_code')%>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE))

biplot_scores <- as.data.frame(vegan::scores(cca_ord, display = "bp")) %>%
  rownames_to_column("Variable")
biplot_scores$Variable <- recode(
  biplot_scores$Variable,
  "subset_meta$time_pointAntibiotics" = "Antibiotics",
  "subset_meta$PopulationWarm" = "Warm population",
  "subset_meta$time_pointAntibiotics:subset_meta$PopulationWarm" =
    "Interaction Wam population"
)

CAP_df %>%
  group_by(type, time_point) %>%
  mutate(x_cen = mean(CAP1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(CAP2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(., aes(x=CAP1,y=CAP2, color=Population,shape = time_point)) +
  scale_color_manual(values=c('#008080', "#d57d2c")) +
  scale_fill_manual(values=c('#00808050', "#d57d2c50")) +
  geom_point(size=2) +
  geom_hline(yintercept = 0, linetype = "dotted") + 
  geom_vline(xintercept = 0, linetype = "dotted") +
  geom_segment(aes(x=x_cen, y=y_cen, xend=CAP1, yend=CAP2), alpha=0.2) +
  geom_segment(data = biplot_scores, aes(x = 0, y = 0, xend = CAP1, yend = CAP2),
               inherit.aes = FALSE, arrow = arrow(length = unit(0.2, "cm")), color = "black") +
  geom_text(data = biplot_scores, aes(x = CAP1, y = CAP2, label = Variable),
            inherit.aes = FALSE, color = "black", vjust = -0.5, hjust = 0.5)+
  theme_classic()

6.3 Differences in MAGs and functional capacity

6.3.1 Cold population

6.3.1.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Cold" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Cold" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

subset_meta <- sample_metadata %>%
  filter(Population == "Cold" & time_point %in% c("Antibiotics","Acclimation"))

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",subset_meta$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% 
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 9 × 2
# Rowwise: 
  phylum               sample_count
  <chr>                       <int>
1 p__Bacillota_A                 43
2 p__Bacteroidota                15
3 p__Bacillota                    8
4 p__Pseudomonadota               8
5 p__Cyanobacteriota              3
6 p__Verrucomicrobiota            2
7 p__Bacillota_B                  1
8 p__Bacillota_C                  1
9 p__Fusobacteriota               1

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)

# A tibble: 82 × 7
# Rowwise: 
   genome                phylum               class               order                 family                 genus          species
   <chr>                 <chr>                <chr>               <chr>                 <chr>                  <chr>          <chr>  
 1 AH1_2nd_1:bin_000006  p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Kineothrix  s__    
 2 AH1_2nd_19:bin_000005 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Roseburia   s__    
 3 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia s__    
 4 AH1_2nd_2:bin_000003  p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__JAAYNV01    s__    
 5 AH1_2nd_10:bin_000007 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__                    g__            s__    
 6 AH1_2nd_15:bin_000043 p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Bacteroidaceae      g__Bacteroides s__    
 7 AH1_2nd_10:bin_000041 p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Rikenellaceae       g__Alistipes   s__    
 8 AH1_2nd_8:bin_000007  p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma      s__    
 9 AH1_2nd_14:bin_000015 p__Bacillota_A       c__Clostridia       o__Christensenellales f__UBA3700             g__            s__    
10 AH1_2nd_5:bin_000001  p__Bacillota_C       c__Negativicutes    o__Selenomonadales    f__                    g__            s__    
# ℹ 72 more rows

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 4 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus         species            
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>         <chr>              
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus    s__Proteus cibarius
2 LI1_2nd_7:bin_000001  p__Bacillota_A    c__Clostridia          o__Clostridiales      f__Clostridiaceae     g__Sarcina    s__                
3 AH1_2nd_7:bin_000055  p__Bacillota      c__Bacilli             o__Mycoplasmatales    f__Mycoplasmoidaceae  g__Ureaplasma s__                
4 AH1_2nd_13:bin_000025 p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA3700            g__           s__                

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>

6.3.1.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Cold" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output_antib_accl_2803 = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output_antib_accl_2803, 
     file = "data/ancombc_antib_accl_2803.Rdata")

load("data/ancombc_antib_accl_2803.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output_antib_accl_2803$res %>%
  dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
  filter(p_time_pointAntibiotics < 0.05) %>%
  dplyr::arrange(p_time_pointAntibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_pointAntibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) + 
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

MAGs enriched in acclimation on the left side and in antibiotics on the right side

Phyla of the significant MAGs

Acclimation

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

            phylum sample_count
1     Bacteroidota           14
2      Bacillota_A            5
3        Bacillota            1
4 Campylobacterota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  select(phylum, genus, species)

             phylum           genus species
1      Bacteroidota     Bacteroides        
2  Campylobacterota  Helicobacter_J        
3       Bacillota_A                        
4      Bacteroidota     Odoribacter        
5         Bacillota                        
6      Bacteroidota     Bacteroides        
7      Bacteroidota     Bacteroides        
8      Bacteroidota     Phocaeicola        
9      Bacteroidota     Bacteroides        
10     Bacteroidota     Odoribacter        
11     Bacteroidota     Phocaeicola        
12     Bacteroidota     Bacteroides        
13     Bacteroidota     Bacteroides        
14     Bacteroidota Parabacteroides        
15     Bacteroidota                        
16     Bacteroidota Parabacteroides        
17      Bacillota_A                        
18      Bacillota_A    Lacrimispora        
19     Bacteroidota       Alistipes        
20      Bacillota_A                        
21      Bacillota_A                        

Antibiotics

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

             phylum sample_count
1       Bacillota_A            2
2         Bacillota            1
3    Pseudomonadota            1
4 Verrucomicrobiota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  select(phylum, genus, species)

             phylum          genus species
1       Bacillota_A     MGBC116941        
2         Bacillota        CAG-345        
3       Bacillota_A Intestinimonas        
4 Verrucomicrobiota                       
5    Pseudomonadota                       

6.3.1.3 Differences in the functional capacity

sample_sub <- sample_metadata %>%
  filter(Population == "Cold" & time_point %in% c("Acclimation", "Antibiotics"))
genome_counts_filtered <- genome_counts_filt %>% 
    select(one_of(c("genome",sample_sub$sample))) 
  
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

genome_counts_filtered_row <- genome_counts_filtered %>%
  mutate_at(vars(-genome),~./sum(.)) %>% 
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)

6.3.1.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.335 0.0323
2 Antibiotics 0.247 0.136 

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.94311, p-value = 0.09178

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(8,9)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(time_point) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
  scale_fill_manual(values=c("#003a45",'#008080')) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="time_point")

6.3.1.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.329 0.0319
2 Antibiotics 0.256 0.126 

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.95736, p-value = 0.2323

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(8,9)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	time_point
B06	0.68059020	0.47329940	Organic anion biosynthesis	0.20729080	Acclimation
B02	0.59930820	0.41576970	Amino acid biosynthesis	0.18353850	Acclimation
D02	0.38530780	0.20616790	Polysaccharide degradation	0.17913990	Acclimation
S03	0.27145162	0.09403129	Spore	0.17742033	Acclimation
B01	0.84159250	0.69078390	Nucleic acid biosynthesis	0.15080860	Acclimation
B07	0.44505530	0.30423320	Vitamin biosynthesis	0.14082210	Acclimation
B08	0.44618700	0.32094170	Aromatic compound biosynthesis	0.12524530	Acclimation
D09	0.25519710	0.13446900	Antibiotic degradation	0.12072810	Acclimation
B04	0.54481600	0.42941310	SCFA biosynthesis	0.11540290	Acclimation
D03	0.29173710	0.20687250	Sugar degradation	0.08486460	Acclimation
D05	0.28853350	0.22270070	Amino acid degradation	0.06583280	Acclimation
B10	0.05960097	0.03635236	Antibiotic biosynthesis	0.02324861	Acclimation

6.3.2 Warm population

6.3.2.1 Structural zeros

acclimation_samples <- sample_metadata %>%
  filter(Population == "Warm" & time_point == "Acclimation") %>% 
  dplyr::select(sample) %>%
  pull()

antibiotics_samples <- sample_metadata %>%
  filter(Population == "Warm" & time_point == "Antibiotics") %>% 
  dplyr::select(sample) %>% pull()

sample_sub <- sample_metadata %>%
  filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics"))

structural_zeros <- genome_counts_filt %>% 
  select(one_of(c("genome",sample_sub$sample))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>% 
   rowwise() %>% 
   mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>% 
   mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>% 
   mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>% 
   mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>% 
   filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE)  %>% 
   mutate(present = case_when(
      all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
      !all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
      !all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
      TRUE ~ NA_character_
    )) %>%
   mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
   dplyr::select(genome, present, average) %>%
   left_join(genome_metadata, by=join_by(genome==genome)) %>%
   arrange(present,-average)

Acclimation

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count)) 

# A tibble: 12 × 2
# Rowwise: 
   phylum               sample_count
   <chr>                       <int>
 1 p__Bacillota_A                 64
 2 p__Bacteroidota                32
 3 p__Bacillota                   16
 4 p__Cyanobacteriota              9
 5 p__Pseudomonadota               8
 6 p__Bacillota_B                  3
 7 p__Desulfobacterota             3
 8 p__Bacillota_C                  2
 9 p__Verrucomicrobiota            2
10 p__Actinomycetota               1
11 p__Campylobacterota             1
12 p__Elusimicrobiota              1

structural_zeros %>% 
  filter(present=="acclimation") %>% 
  select(1,5:10)

# A tibble: 142 × 7
# Rowwise: 
   genome                phylum               class               order                 family                 genus            species
   <chr>                 <chr>                <chr>               <chr>                 <chr>                  <chr>            <chr>  
 1 AH1_2nd_11:bin_000008 p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia   s__    
 2 LI1_2nd_8:bin_000013  p__Verrucomicrobiota c__Verrucomicrobiae o__Verrucomicrobiales f__Akkermansiaceae     g__Akkermansia   s__    
 3 AH1_2nd_19:bin_000005 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Roseburia     s__    
 4 AH1_2nd_20:bin_000009 p__Bacillota_A       c__Clostridia       o__Lachnospirales     f__Lachnospiraceae     g__Kineothrix    s__    
 5 AH1_2nd_18:bin_000011 p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Coprobacillaceae    g__Coprobacillus s__    
 6 LI1_2nd_2:bin_000001  p__Bacillota_A       c__Clostridia       o__Oscillospirales    f__Ruminococcaceae     g__              s__    
 7 AH1_2nd_6:bin_000035  p__Bacillota         c__Bacilli          o__RF39               f__UBA660              g__Faecisoma     s__    
 8 LI1_2nd_4:bin_000003  p__Bacillota         c__Bacilli          o__Erysipelotrichales f__Erysipelotrichaceae g__Dielma        s__    
 9 AH1_2nd_1:bin_000039  p__Bacteroidota      c__Bacteroidia      o__Bacteroidales      f__Marinifilaceae      g__Odoribacter   s__    
10 AH1_2nd_2:bin_000018  p__Elusimicrobiota   c__Elusimicrobia    o__Elusimicrobiales   f__Elusimicrobiaceae   g__UBA1436       s__    
# ℹ 132 more rows

Antibiotics

structural_zeros %>% 
  filter(present=="antibiotics") %>% 
  select(1,5:10)

# A tibble: 7 × 7
# Rowwise: 
  genome                phylum            class                  order                 family                genus              species                       
  <chr>                 <chr>             <chr>                  <chr>                 <chr>                 <chr>              <chr>                         
1 AH1_2nd_7:bin_000003  p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales   f__Enterobacteriaceae g__Proteus         s__Proteus cibarius           
2 LI1_2nd_8:bin_000030  p__Actinomycetota c__Actinomycetia       o__Mycobacteriales    f__Mycobacteriaceae   g__Corynebacterium s__                           
3 LI1_2nd_8:bin_000077  p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__Parabacteroides sp003480915
4 LI1_2nd_5:bin_000032  p__Bacillota_A    c__Clostridia          o__Christensenellales f__UBA1242            g__Caccovivens     s__                           
5 AH1_2nd_18:bin_000013 p__Bacteroidota   c__Bacteroidia         o__Bacteroidales      f__Tannerellaceae     g__Parabacteroides s__                           
6 LI1_2nd_5:bin_000069  p__Bacillota      c__Bacilli             o__RF39               f__UBA660             g__                s__                           
7 AH1_2nd_20:bin_000061 p__Bacillota      c__Bacilli             o__Lactobacillales    f__Enterococcaceae    g__Enterococcus    s__                           

Community elements differences in structural zeros

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)

# A tibble: 20 × 3
   trait p_value p_adjust
   <chr>   <dbl>    <dbl>
 1 B0220 0.00618   0.0497
 2 B0708 0.00280   0.0497
 3 B0709 0.00662   0.0497
 4 B1012 0.00618   0.0497
 5 D0203 0.00559   0.0497
 6 D0205 0.00280   0.0497
 7 D0206 0.00280   0.0497
 8 D0207 0.00280   0.0497
 9 D0210 0.00280   0.0497
10 D0211 0.00662   0.0497
11 D0213 0.00685   0.0497
12 D0305 0.00685   0.0497
13 D0308 0.00685   0.0497
14 D0610 0.00618   0.0497
15 D0706 0.00662   0.0497
16 D0902 0.00618   0.0497
17 D0905 0.00618   0.0497
18 D0908 0.00280   0.0497
19 D0911 0.00618   0.0497
20 S0104 0.00685   0.0497

6.3.2.2 Differential abundance analysis: composition

phylo_samples <- sample_metadata %>%
  filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics") )%>% 
  column_to_rownames("sample") %>% 
  sample_data()
phylo_genome <- genome_counts_filt %>% 
  filter(!genome %in% structural_zeros$genome) %>% 
  select(one_of(c("genome",rownames(phylo_samples)))) %>%
  filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
  select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
  column_to_rownames("genome") %>% 
  otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>% 
  filter(genome %in% rownames(phylo_genome)) %>% 
  column_to_rownames("genome") %>% 
  dplyr::select(domain,phylum,class,order,family,genus,species) %>% 
  as.matrix() %>% 
  tax_table()

physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)

ancom_rand_output_2803 = ancombc2(data = physeq_genome_filtered, 
                  assay_name = "counts", 
                  tax_level = NULL,
                  fix_formula = "time_point",
                  rand_formula = "(1|individual)",
                  p_adj_method = "holm", 
                  pseudo_sens = TRUE,
                  prv_cut =0, 
                  lib_cut = 0, 
                  s0_perc = 0.05,
                  group = NULL, 
                  struc_zero = FALSE, 
                  neg_lb = FALSE,
                  alpha = 0.05, 
                  n_cl = 2, 
                  verbose = TRUE,
                  global = FALSE, 
                  pairwise = FALSE, 
                  dunnet = FALSE, 
                  trend = FALSE,
                  iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
                  em_control = list(tol = 1e-5, max_iter = 100),
                  lme_control = lme4::lmerControl(),
                  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100), 
                  trend_control = NULL)

save(ancom_rand_output_2803, 
     file = "data/ancombc_antib_accl_warm_2803.Rdata")

load("data/ancombc_antib_accl_warm_2803.Rdata")

taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
  rownames_to_column(., "taxon") %>%
  mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))

ancombc_rand_table_mag <- ancom_rand_output_2803$res %>%
  dplyr::select(taxon, lfc_time_pointAntibiotics, p_time_pointAntibiotics) %>%
  filter(p_time_pointAntibiotics < 0.05) %>%
  dplyr::arrange(p_time_pointAntibiotics) %>%
  merge(., taxonomy, by="taxon") %>%
  mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
  dplyr::arrange(lfc_time_pointAntibiotics)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", ""))  %>%
  right_join(taxonomy, by=join_by(phylum == phylum)) %>%
  dplyr::select(phylum, colors) %>%
  mutate(colors = str_c(colors, "80"))  %>% #add 80% alpha
    unique() %>%
    dplyr::arrange(phylum)

tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
  
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
    dplyr::arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

ancombc_rand_table_mag%>%
      mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_pointAntibiotics, y=forcats::fct_reorder(genome,lfc_time_pointAntibiotics), fill=phylum)) + 
  geom_col() + 
  scale_fill_manual(values=tax_color) + 
  geom_hline(yintercept=0) + 
  theme(panel.background = element_blank(),
        axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 8),
        axis.title = element_text(size = 14, face = "bold"),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14, face = "bold"),
        legend.position = "right", legend.box = "vertical")+
  xlab("log2FoldChange") + 
  ylab("Species")+
  guides(fill=guide_legend(title="Phylum"))

MAGs enriched in acclimation on the left side and in antibiotics on the right side

Phyla of the significant MAGs Acclimation

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics<0)  %>% 
  select(phylum, genus, species)

        phylum       genus species
1 Bacteroidota Phocaeicola        

Antibiotics

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  count(phylum, name = "sample_count") %>%
  arrange(desc(sample_count))  

          phylum sample_count
1   Bacteroidota            2
2      Bacillota            1
3    Bacillota_A            1
4    Chlamydiota            1
5 Pseudomonadota            1

ancombc_rand_table_mag%>%
  filter(lfc_time_pointAntibiotics>0)  %>% 
  select(phylum, genus, species)

          phylum       genus species
1   Bacteroidota Bacteroides        
2    Bacillota_A  MGBC116941        
3   Bacteroidota Bacteroides        
4      Bacillota  Ureaplasma        
5 Pseudomonadota                    
6    Chlamydiota                    

6.3.2.3 Differences in the functional capacity

sample_sub <- sample_metadata %>%
  filter(Population == "Warm" & time_point %in% c("Acclimation", "Antibiotics"))

genome_counts_filtered <- genome_counts_filt %>% 
    select(one_of(c("genome",sample_sub$sample)))

GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filtered$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>% 
  select_if(~ !is.numeric(.) || sum(.) != 0)

GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)

GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)

genome_counts_filtered_row <- genome_counts_filtered %>%
  mutate_at(vars(-genome),~./sum(.))  %>%
  column_to_rownames(., "genome") 

GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_filtered_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_filtered_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_filtered_row,GIFT_db)

6.3.2.3.1 Community elements differences:

MCI

GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.355 0.0228
2 Antibiotics 0.253 0.0674

MCI <- GIFTs_elements_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.88114, p-value = 0.03317

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 70, p-value = 0.0003291
alternative hypothesis: true location shift is not equal to 0

element_gift <- GIFTs_elements_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  left_join(., sample_metadata[c(8,9)], by=join_by("sample"=="sample"))

uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)

significant_elements <- element_gift %>%
    pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
    group_by(trait) %>%
    summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
    mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
    filter(p_adjust < 0.05)%>%
  left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))

element_gift_t <- element_gift  %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "trait")

element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>% 
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric)  %>%
  rownames_to_column(., "sample")%>% 
  left_join(., sample_metadata[c(8,9)], by = join_by(sample == sample))

difference_table <- element_gift_filt %>%
  select(-sample) %>%
  group_by(time_point) %>%
  summarise(across(everything(), mean)) %>%
  t() %>%
  row_to_names(row_number = 1) %>%
  as.data.frame() %>%
  mutate_if(is.character, as.numeric) %>%
  rownames_to_column(., "Elements") %>%
  left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>% 
  arrange(Function) %>% 
  mutate(Difference=Acclimation-Antibiotics)%>% 
  mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) 

difference_table %>%
  ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) + 
  geom_col() +
  scale_fill_manual(values=c("#003a45",'#008080')) +
  geom_hline(yintercept=0) + 
  coord_flip()+
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        legend.position = "right", 
        panel.background = element_blank(),
          panel.grid.major = element_line(size = 0.15, linetype = 'solid',
                                colour = "grey"))+
  xlab("Microbial Functional Capacity") + 
  ylab("Mean difference")+
  labs(fill="time_point")

6.3.2.3.2 Community functions differences

MCI

GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(time_point) %>%
  summarise(MCI = mean(value), sd = sd(value))

# A tibble: 2 × 3
  time_point    MCI     sd
  <chr>       <dbl>  <dbl>
1 Acclimation 0.348 0.0158
2 Antibiotics 0.268 0.0586

MCI <- GIFTs_functions_community %>%
  rowMeans() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) 

shapiro.test(MCI$value)

    Shapiro-Wilk normality test

data:  MCI$value
W = 0.85791, p-value = 0.01419

wilcox.test(value ~ time_point, data=MCI)

    Wilcoxon rank sum exact test

data:  value by time_point
W = 69, p-value = 0.0005759
alternative hypothesis: true location shift is not equal to 0

function_gift <- GIFTs_functions_community %>% 
  as.data.frame() %>% 
  rownames_to_column(., "sample") %>% 
  merge(., sample_metadata[c(8,9)], by="sample")

unique_funct_db<- GIFT_db[c(3,4,5)] %>% 
  distinct(Code_function, .keep_all = TRUE)

significant_functional <- function_gift %>%
  pivot_longer(-c(sample,time_point), names_to = "trait", values_to = "value") %>%
  group_by(trait) %>%
  summarise(p_value = wilcox.test(value ~ time_point)$p.value) %>%
  mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
  filter(p_adjust < 0.05)%>%
  left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))

Code_function	Acclimation	Antibiotics	Function	Difference	time_point
D02	0.42716670	0.21793910	Polysaccharide degradation	0.20922760	Acclimation
B02	0.62784110	0.47187490	Amino acid biosynthesis	0.15596620	Acclimation
D03	0.34000130	0.18813800	Sugar degradation	0.15186330	Acclimation
B08	0.44482920	0.33165890	Aromatic compound biosynthesis	0.11317030	Acclimation
D05	0.33300830	0.22491860	Amino acid degradation	0.10808970	Acclimation
D09	0.24738110	0.14240830	Antibiotic degradation	0.10497280	Acclimation
B06	0.70894950	0.61231610	Organic anion biosynthesis	0.09663340	Acclimation
B09	0.02947309	0.01568554	Metallophore biosynthesis	0.01378755	Acclimation